home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / sqrt.c < prev    next >
C/C++ Source or Header  |  1990-11-23  |  1KB  |  48 lines
  1. /* SQRT function for Sozobon C.                        */
  2. /* Copyright © David Brooks, 1989 All Rights Reserved            */
  3. /*                                    */
  4. /* The method is from Hart et al: "Computer Approximations".  We reduce    */
  5. /* the range to [0.25..1.0], then apply the polynomial:            */
  6. /* 0.25927688 + 1.0520212n - 0.31632214n^2, giving 2.3 digits         */
  7. /* approximation.  Then we do two Newton's iterations (giving ~9 digits    */
  8. /* precision).                                */
  9. /*                                    */
  10. /* A negative argument gives a return value of 0.0 and sets EDOM.    */
  11.  
  12. #include <fplib.h>
  13. #include <errno.h>
  14.  
  15. /* Sadly we can't use register for objects of type fstruct. */
  16.  
  17. float sqrt(n)
  18. fstruct n;
  19. {    register int exp;        /* Separate unbiased exponent.    */
  20.     fstruct s;
  21.  
  22.     if (n.sc[3] == 0)        /* Check for 0.0        */
  23.         return 0.0;
  24.  
  25.     if ((exp = n.sc[3]) < 0)    /* Test argument's sign.    */
  26.     {    errno = EDOM;
  27.         return 0.0;
  28.     }
  29.  
  30.     exp -= BIAS;            /* Get unbiased exponent.    */
  31.     if (exp & 1)
  32.     {    n.sc[3] = BIAS-1;    /* Convert to 0.25..0.5.    */
  33.         ++exp;
  34.     }
  35.     else
  36.         n.sc[3] = BIAS;        /* Or convert to 0.5..1.    */
  37.  
  38.     s.f = 0.25927688 + n.f * (1.0520212 + n.f * -0.31632214);
  39.  
  40.     s.f = s.f + n.f / s.f;        /* One Newton.            */
  41.     s.sc[3] -= 1;            /* (here's the divide by 2).    */
  42.     s.f = s.f + n.f / s.f;        /* The other Newton.        */
  43.     s.sc[3] += ((unsigned int)exp >> 1) - 1;
  44.                     /* Divide by 2 and insert
  45.                        half the original exponent.    */
  46.     return s.f;
  47. }
  48.